
merged = read.csv(file.choose())

merged[merged$period=='current',][5:6]= 'Current'

# Area change analysis
Area = data.frame(
  
  # retain all needed GCMs
  GCM = c('Current',
          'HadGEM3','HadGEM3',
          'ACCESS','ACCESS',
          'MIROC6','MIROC6',
          'CMCC','CMCC',
          'INM','INM',
          'Mean','Mean'),
  
  period=c('Current',
           '2041-2060','2061-2080',
           '2041-2060','2061-2080',
           '2041-2060','2061-2080',
           '2041-2060','2061-2080',
           '2041-2060','2061-2080',
           '2041-2060','2061-2080'),
  
  # follow same order as in GCM colum created above
  suitable_area=c(nrow(filter(merged,GCM=='Current'& binary=='Suitable')),
                  nrow(filter(merged,GCM=='HadGEM3'& period=='2041-2060' & binary =='Suitable')),
                  nrow(filter(merged,GCM=='HadGEM3'& period=='2061-2080' & binary =='Suitable')),
                  nrow(filter(merged,GCM=='ACCESS'& period=='2041-2060' & binary =='Suitable')),
                  nrow(filter(merged,GCM=='ACCESS'& period=='2061-2080' & binary =='Suitable')),
                  nrow(filter(merged,GCM=='MIROC6'& period=='2041-2060' & binary =='Suitable')),
                  nrow(filter(merged,GCM=='MIROC6'& period=='2061-2080' & binary =='Suitable')),
                  nrow(filter(merged,GCM=='CMCC'& period=='2041-2060' & binary =='Suitable')),
                  nrow(filter(merged,GCM=='CMCC'& period=='2061-2080' & binary =='Suitable')),
                  nrow(filter(merged,GCM=='INM'& period=='2041-2060' & binary =='Suitable')),
                  nrow(filter(merged,GCM=='INM'& period=='2061-2080' & binary =='Suitable')),
                  nrow(filter(merged,GCM=='Mean'& period=='2041-2060' & binary =='Suitable')),
                  nrow(filter(merged,GCM=='Mean'& period=='2061-2080' & binary =='Suitable'))
  ))

# percantage change in suitability area relative to the current suitability area size
Area = Area %>% mutate(prop=(suitable_area-Area[1,3])/Area[1,3])

Area$period = factor(Area$period,levels = c('Current','2041-2060','2061-2080'))

Area$GCM = factor(Area$GCM,levels=c('Current','HadGEM3','ACCESS','MIROC6','CMCC','INM','Mean'))

Area$scenario = 'SSP2-4.5'
#-------------------------------------------------------------------------------
merged1 = read.csv(file.choose())

merged1[merged1$period=='current',][5:6]= 'Current'

# Area change analysis
Area1 = data.frame(
  
  # retain all needed GCMs
  GCM = c('Current',
          'HadGEM3','HadGEM3',
          'ACCESS','ACCESS',
          'MIROC6','MIROC6',
          'CMCC','CMCC',
          'INM','INM',
          'Mean','Mean'),
  
  period=c('Current',
           '2041-2060','2061-2080',
           '2041-2060','2061-2080',
           '2041-2060','2061-2080',
           '2041-2060','2061-2080',
           '2041-2060','2061-2080',
           '2041-2060','2061-2080'),
  
  # follow same order as in GCM colum created above
  suitable_area=c(nrow(filter(merged1,GCM=='Current'& binary=='Suitable')),
                  nrow(filter(merged1,GCM=='HadGEM3'& period=='2041-2060' & binary =='Suitable')),
                  nrow(filter(merged1,GCM=='HadGEM3'& period=='2061-2080' & binary =='Suitable')),
                  nrow(filter(merged1,GCM=='ACCESS'& period=='2041-2060' & binary =='Suitable')),
                  nrow(filter(merged1,GCM=='ACCESS'& period=='2061-2080' & binary =='Suitable')),
                  nrow(filter(merged1,GCM=='MIROC6'& period=='2041-2060' & binary =='Suitable')),
                  nrow(filter(merged1,GCM=='MIROC6'& period=='2061-2080' & binary =='Suitable')),
                  nrow(filter(merged1,GCM=='CMCC'& period=='2041-2060' & binary =='Suitable')),
                  nrow(filter(merged1,GCM=='CMCC'& period=='2061-2080' & binary =='Suitable')),
                  nrow(filter(merged1,GCM=='INM'& period=='2041-2060' & binary =='Suitable')),
                  nrow(filter(merged1,GCM=='INM'& period=='2061-2080' & binary =='Suitable')),
                  nrow(filter(merged1,GCM=='Mean'& period=='2041-2060' & binary =='Suitable')),
                  nrow(filter(merged1,GCM=='Mean'& period=='2061-2080' & binary =='Suitable'))
  ))

# percantage change in suitability area relative to the current suitability area size
Area1 = Area1 %>% mutate(prop=(suitable_area-Area1[1,3])/Area1[1,3])

Area1$period = factor(Area$period,levels = c('Current','2041-2060','2061-2080'))

Area1$GCM = factor(Area1$GCM,levels=c('Current','HadGEM3','ACCESS','MIROC6','CMCC','INM','Mean'))

Area1$scenario = 'SSP5-8.5'

# merge data
Area2 = rbind(Area,Area1)

Area2$period = ifelse(Area2$period=='2041-2060','2050s',ifelse(Area2$period=='2061-2080',
                                                               '2070s',Area2$period))
Area2[Area2$period=='Current'][1]= 'Current'

# extract current data only
current_1 = filter(Area,GCM=='Current') %>% dplyr::select(-c('GCM'))


# plotting the results
p2 = ggplot(data=filter(Area2,GCM!='Current'),aes(period,suitable_area,fill=scenario))+
    geom_col(width=0.4,position = 'dodge')+
    geom_col(data=current_1,aes(period,suitable_area),fill = 'green3',width=0.4)+
    geom_text(aes(label=suitable_area),position=position_dodge(width=0.5),vjust=-0.5,size=2)+
    geom_text(data=current_1,aes(label=suitable_area),vjust=-0.5,size=2)+
    labs(y=expression(Total~potential~suitable~area~'('~km^2~')'),x='',fill='')+
    theme_bw()+
    theme(legend.position = 'none')+
    scale_fill_manual(values= c('indianred1','cyan3'))+
    scale_y_continuous(limits = c(0,4000))+
    facet_wrap(~GCM)

#-------------------------------------------------------------------------------
p3 = ggplot(filter(Area2,GCM!='Current'),aes(period,prop,fill=scenario))+
  geom_col(width=0.4,position = 'dodge')+
  geom_text(aes(label=scales::percent(prop,accuracy = 0.1)),position=position_dodge(width=0.5),vjust=1.2,size=2)+
  scale_y_continuous(labels = scales::percent,limits = c(-1.1,0))+
  labs(y=expression(Change~'in'~area~suitability),x='',fill='')+
  theme_bw()+
  #scale_y_continuous(limits = c(-1,0.4))+ # only applies for Arabica
  theme(legend.position = 'bottom')+
  facet_wrap(~GCM)

p4 = p2/p3

p4

ggsave('Area_change_Robusta_analysis_ssp245_ssp585.png',p4,height=6,width=9,units='in',dpi=320)

rm(list=ls())
dev.off()
